Thermodynamics and Kinetics of a Go Protein-Like 
Heteropolymer Model with Two-State Folding Characteristics 
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We present results of Monte Carlo computer simulations of a coarse-grained hydrophobic-polar 
Go-like heteropolymer model and discuss thermodynamic properties and kinetics of an exemplified 
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I. INTRODUCTION 

Spontaneous protein folding is a dynamic process, 
which starts after the generation of the DNA-encoded 
amino acid sequence in the ribosome and is in many cases 
finished, when the functional conformation, the native 
fold, is formed. As this process takes microseconds to 
seconds, a dynamical computational analysis of an ap- 
propriate microscopic model, which could lead to a bet- 
ter understanding of the conformational transitions ac- 
companying folding [1], is extremely demanding. Since 
protein folding is a thermodynamic process at finite tem- 
perature, a certain folding trajectory in the free-energy 
landscape is influenced by Brownian collisions with sur- 
rounding solvent molecules. Therefore, it is more favor- 
able to study the kinetics of the folding process by aver- 
aging over an appropriate ensemble of trajectories. 

A significant problem is that the complexity of de- 
tailed semiclassical microscopic models based on force- 
fields and solvent parameter sets (or explicit solvent) 
rules out molecular dynamics (MD) in many cases and, 
therefore, Markovian Monte Carlo (MC) dynamics is a 
frequently used method for such kinetic studies. It is ob- 
vious, however, that the time scale provided by MC is 
not directly comparable with the time scale of the fold- 
ing process. It is widely believed that the folding path 
a protcir 
1, i.e., _ 

between residues and, therefore, long-range correlations 
and memory effects can significantly influence the kinet- 
ics. 

A few years ago, experimental evidence was found that 
classes of proteins show particular simple folding charac- 
teristics, single exponential and two-state folding [1, 0]. 



of a protein is strongly correlated with contact order- 
ing Ql, i.e., the order of the successive contact formation 
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In the two-state folding process, which is in the focus of 
the present study, the peptide is either in an unfolded, 
denatured state or it possesses a native-like, folded struc- 
ture. In contrast to the barrier-free single-exponential 
folding, there exists an unstable transition state to be 
passed in the two-state folding process. Due to the com- 
paratively simple folding characteristics, strongly simpli- 
fied, effective models were established. Knowledge-based 
models of Go-like type la, [ 
mcrous recent studies 
In Go-like models the native fold must be known and is 
taken as input for the energy function. The energy of 
an actual conformation depends on its structural devia- 
tion from the native fold (e.g., by counting the number of 
already established native contacts). By definition, the 
energy is minimal, if conformation and native fold are 
identical in all degrees of freedom involved in the model. 
The simplicity of the model entails reduced computa- 
tional complexity and also MD simulations, e.g., based on 
Langevin dynamics [ifj, can successfully be performed. 



In this paper, we follow a different approach. We also 
study a Go-like model, but it is based on a minimalistic 
coarse-grained hydrophobic-polar representation of the 
heteropolymer. The basic idea behind coarse-graining is 
the introduction of a mesoscopic length scale, i.e., the 
reduction of microscopic details, in order to classify het- 
cropolymers with respect to their folding characteristics. 
The assumption is that, if there is some sort of universal- 
ity in folding processes, then an effective model should 
allow for a general description of the qualitative fold- 
ing behavior. In fact, in a recent work we could show 
that two-state folding, folding through intermediates, 
and metastability are inherent tertiary fol ding processes 
which can also be found in folding studies [ijj, [IS] of the 
simple hydrophobic-polar AB model [2l| . In this con- 
text, it is quite interesting that also secondary structures 
are intrinsic geometries of pol yme r-like objects, even on 
a mesoscopic scale [H, HH EE HI] ■ 
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II. MODEL AND METHODS 

In the following, we consider the hydrophobic- 
polar heteropolymer sequence S = A 3 B 2 AB 2 ABA- 
B 2 AB ABABA, where the A's indicate hydrophobic 
monomers and the B's polar (hydrophilic) residues. For 
our comparative model study, we employ the physically 
motivated AB model l2ll l22l ] and a knowledge-based 
model of Go type [H, Irl 111] , which is referred to as the 
GoL (Go-like) model throughout the paper. 

In the AB model, the energy of a conformation R = 
{ri, ... tat} with N monomers, where Ti denotes the spa- 
tial location of the ith monomer, is given by 



N-2 




cos fil 



C(<Ji,<jj) 



(1) 



where the first term is the bending energy and the sum 
runs over the (N — 2) bending angles dk (0 < fit < tt) 
between monomers fc, k + 1, and k + 2 of successive cova- 
lent bond vectors which have length unity. The Lennard- 
Jones type potentials depend on the types of the inter- 
acting monomers (ai = A, B) and on their spatial dis- 
tance Tij and represent the influence of the specific AB 
sequence on the energy. The long-range behavior is at- 
tractive for pairs of like monomers and repulsive for AB 
pairs of monomers: 




a l ,<jj = A, 
a tl aj = B, 
<j l ^ a j . 



(2) 



In this model, the AB heteropolymer with sequence S ex- 
periences a hydrophobic collapse transition which is sig- 
nalized by the peak in the specific-heat curve plotted in 
Fig. [T] The transition between random coils and native- 
like hydrophobic-core conformations is a (pseudo) phase 
separation process and the folding transition of two-state 
(folded/unfolded) type. In particular, for such systems 
it is known from model studies of realistic amino acid 
sequences that knowledge-based Go type models reveal 
kinetic as pec ts of folding and unfolding processes reason- 
ably well [H EU ■ 

For performing kinetic studies of S in the simpli- 
fied hydrophobic-polar approach as well, we use the 
(putative) global-energy minimum, identified in energy- 
landscape paving (ELP) minimizations [27| of the AB 
model [28] . as input for the definition of a hydrophobic- 
polar GoL model. The (putative) native conformation 
n/ -* is shown in Fig. [5] and its energy is Eab ~ —19.3 
in the units of the AB model JT]). In GoL models, the 
"energy" of a given conformation is related to its similar- 
ity with the ground state. This means, "energy" in the 
GoL model plays rather the role of a similarity or "or- 
der" parameter and is, therefore, not a potential energy 
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FIG. 1: Specific heats for t he p eptide with sequence S as 
obtained with the AB model [2a ] and the gauged GoL model 
as defined in Eq. @. 




FIG. 2: Global-energy minimum of sequence S in the AB 
model with -Eab ~ —19.3. Dark monomers are hydrophobic 
(A) and light residues polar (B). 



in the usual physical sense (as there is no physical force 
associated with it). Denoting the (N — 2) bending angles 
of the global-energy minimum conformation by #1 , its 
( N — 3) torsional angles as (p^ , and the inter-monomer 
distances by , we define the GoL model according to 
the representation in Ref. [ll| as: 
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The last two sums run over all pairs of nonbonded 
monomers. In the case the pair forms a native 

contact, i.e., ry < r cut and < r cut , the monomers 
experience a short-range 10-12 Lennard-Jones attrac- 
tion, while for nonnative contacts an overall repulsive 
1/r 12 contribution is taken into account. The constants 
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FIG. 3: Fluctuations of the similarity parameter g and the 
angular overlap parameter Q as functions of the temperature. 
For comparison, the specific heat (dotted line) is also plotted 
into this figure. 

K$ = 20, K^P = 0.5, and = 0.25 weight the rela- 
tive strengths of the angular energy contributions. The 
values were adjusted to have a reasonable coincidence of 
the peak temperature of the specific heat compared to 
the results obtained with the AB model [H[ (see Fig. [1]). 
The free overall energy scale e was chosen such that 

£ AB (R (0) ) = £gol(R (0) ) = -entot, where n tot is the 
total number of native contacts. The definition of a na- 
tive contact requires the introduction of a cutoff radius, 
which we chose as r cut = 1.14. This entails n to t = 20 for 
the native conformation R*- ^, and therefore e w 0.966. 
The choice of the cutoff radius is not very sensitive with 
respect to the results, provided it is not too small. 

For the thermodynamic analyses, we performed par- 
allel tempering simulations [29j, |30j, l3l| and for the ki- 
netic studies, a standard implementation of the Metropo- 
lis Monte Carlo method was used. Chain updates were 
performed employing the spherical-cap updates as de- 
scribed and parametrized in Ref. [2^|. A sweep consists 
of (N — 1) sequential bond vector update trials. 



III. THERMODYNAMICS 

For the study of the folding transition, the introduc- 
tion of an effective parameter is useful which uniquely 
describes the macrostate of the ensemble of heteropoly- 
mer conformations. A widely used measure is the con- 
tact number q(R) which is for a given conformation R 
simply defined as the fraction of the already formed 
native contacts rt(R) in conformation R and the total 
number of native contacts n to t in the final fold, i.e., 
q(R) = 7i(R)/ritot' Then, the statistical ensemble aver- 
age of this quantity (q(R)) at a given temperature char- 
acterizes its macrostate. Roughly, for a two-state folder, 
if (g(R)) > 0.5, native-like conformations are dominat- 
ing the statistical ensemble. If less than half the total 
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FIG. 4: Free-energy landscape F(q) of S for different temper- 
atures close to the transition point. 

number of contacts is formed, the heteropolymer tends 
to reside in the pseudophase of denatured conformations. 
Note that folding transitions are not sharp phase transi- 
tions in the thermodynamic sense as the heteropolymer 
sequence is of finite length and cannot be extended. 

Another suitable parameter is the angular overlap pa- 
rameter [28| which has been proven to be extremely use- 
ful in the characterization of heteropolymer folding chan- 
nels in the AB model [ijj [2(| • If all virtual bond and tor- 
sion angles of two conformations, to be compared with 
each other, coincide, it is unity and between and 1 
otherwise. The advantage is that it is particularly use- 
ful for classifying intermediate or metastable structures 
with stable, but nonnative contacts. In Fig. [3l the fluc- 
tuations of both parameters, i.e., d(q)/dT and d(Q)/dT, 
respectively, are plotted as functions of the temperature 
for the GoL model of sequence S. We clearly see that 
the temperature region of conformational activity as sig- 
nalized by these two "order" parameters coincides with 
the thermally active region indicated by the peak of the 
specific heat, which is also shown for comparison. The 
folding temperature, i.e., the temperature of maximum 
activity, is Tf w 0.36. 

The classification of the heteropolymer with sequence 
S as two-state folder arises from the analysis of the free- 
energy landscape. We assume that q is a suitable pa- 
rameter that describes the macrostate of the system ad- 
equately. Considering this parameter as a constraint, 
we can formally average out the conformational degrees 
of freedom and the probability for a conformation in a 
macrostate with contact parameter q' reads 

P(q') - (%'-g(R))) 

= I JvR6(q' -g(R))e- jEG5L(R)/,CsT , (4) 

where Z denotes the unrestricted partition func- 
tion. The integral measure is simply 2?R = 
ntil^inti^G^+il - !)]■ Expression © can be 
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FIG. 5: Radial distribution functions of S (a) for the lowest- 
energy conformation (see Fig. [2}, (b) at 16 temperatures in 
the interval T £ [0.1, 1.5] (curves from top to bottom) employ- 
ing the AB model, and (c) at 18 temperatures in the interval 
T G [0.1, 1.5] for the GoL model. 



used to define a free energy as function of the q parame- 
ter by 



F(q) = -k B T\np(q). 



(5) 



Since the value of q is a qualitative measure for the 
macrostate the system resides in, the minimum of the 
free energy at a given temperature T is related to the 
dominant macrostate in the canonical ensemble at this 
temperature. Actually, as can be read off from Fig.|H the 
folding transition of the heteropolymer with sequence S 
is a phase-separation process, i.e., at the transition point 
close to Tf « 0.36, the folded and the denatured pseu- 
dophase coexist and the transition state barrier possesses 
a local maximum close to q ~ 0.5, as expected for a typ- 
ical two-state folding characteristics. 



A useful measure of structure formation is the radial 
distribution function and its dependence on temperature. 
Due to translational invariance in the three-dimensional 
space, three degrees of freedom contribute only a volume 
factor to the partition function. Therefore, we utilize this 
by fixing the position of the first monomer, ri = 0. Thus, 
we measure radial distances r of the other monomers with 
respect to the first one and we define the radial distribu- 
tion function as 



G(r) 




(6) 



Note that we have excluded i = 2 as, by definition, the 
virtual covalent bonds are rigid and, therefore, r\2 = 1. 
Actually, from our definition, An J °° drr 2 G{r) = N — 2. 
In Fig. [5l the radial distribution functions of S are shown 
for (a) the lowest-energy conformation and the ensem- 
bles at different temperatures in the (b) AB model and 
(c) GoL model. Although the global-energy minimum 
conformation shown in Fig. [2] does obviously not form 
a regular crystal structure, the Lennard- Jones like in- 
teractions in combination with the rigid virtual bonds 
induce preferable substructures, favoring, e.g., bond an- 
gles of 60°, 90°, and 120° 20]. Therefore, the peaks 
of the radial distances from the first monomer for the 
global energy minimum conformation in Fig.[5fa) can be 
partly explained by these local segments. The first peak 
at r = 2 1/6 « 1.12 is related to the minimum potential 
distance between two nonbonded A monomers (and the 
reference monomer at ri = is of type A). Also the 
location of other peaks can be deduced from Fig. [2] by 
similar geometric arguments. 

The temperature dependence of the peak evolution is 
shown for the AB model in Fig. EJb) and for the GoL 
model in Fig. OUc). As a first result, we see that there is 
nice coincidence not only in the location of the peaks, 
but also in the fact that their is no indication of in- 
termediary, weakly stable conformations. Actually, the 
folding characteristics is similar in both models, at least 
from the thermodynamic point of view. At high temper- 
atures, random-coil conformations dominate: There is 
no significant structuring in the radial distribution func- 
tion. Passing the folding transition temperature, local, 
planar structures form first, before the tertiary, three- 
dimensional ordering towards the native conformation 
occurs. 



IV. KINETICS 

The advantage of the simple GoL model for two-state 
folding is that it enables likewise kinetic studies of folding 
and unfolding events. In fact, this is the main purpose of 
this knowledge-based model, because kinetics studies of 
physically motivated models are typically computation- 
ally extremely demanding. This is unfortunately also 
the case for folding studies employing the AB model. A 
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FIG. 6: Gray-scale coded averages of the probabilities for 
native-contact formation of sequence S as a function of Monte 
Carlo time for the (a) AB and the (b) GoL model, averaged 
over 1 000 folding events. The probability is calculated as a 
combined temporal average over 2 000 MC steps in the en- 
semble of the 1 000 folding events. The darker the bars are, 
the higher is the probability that the associated contact has 
formed. The labels of the pair contacts refer to the numbering 
in the native conformation shown in Fig. m 



striking example is shown in Fig. [5] where for a folding 
event the gray-scale coded average probability of native- 
contact formation is plotted for the original AB model 
[Fig. EJa)] and the GoL model [Fig. E£b)]. The average 
was taken over a "time" window [£mc — At, tyic + At] 



with At — 1 000 MC steps. The darker a bar, the higher 
the probability that the associated contact is formed. As 
the simulation was carried out in the folding regime, the 
probability of native-bond formation increases with the 
number of total MC steps. It is, however, obvious that 
the folding is much slower using the AB model, whereas 
in the GoL model most of the native contacts have al- 
ready been formed after 80000 MC steps. The native 
conformation of the considered sequence possesses partly 
kind of zig-zag structures, or "turns" (see Fig. [2|). The 
folding of these segments is particularly simple and the 
probability of the formation of the associated native con- 
tacts of monomers i with monomers i+2 or i + 3 increases 
much faster than for the other contacts. But even in this 
case, the GoL kinetics is unbeatable. 

This example exhibits the dilemma of physics-based 
models in studying kinetic aspects of structure forma- 
tion by means of computer simulations. It is a notori- 
ously difficult problem, because the time scale of molec- 
ular dynamics is typically too small to see folding events, 
but also Markov chain Monte Carlo dynamics of phys- 
ical models is typically too slow. Just for kinetic as- 
pects, where absolute time scales are widely irrelevant, 
knowledge-based GoL models are an alternative that al- 
low for increasing the sampling efficiency, at least for 
systems with simple folding characteristics (as, for ex- 
ample, two-state folding). Actually, for such models an 
adequate statistical sampling of ensembles close to the 
transition state can be achieved by sophisticated Monte 
Carlo methods where efficiency is typically gained by sim- 
ulating generalized ensembles at the expense of an ar- 
tificial dynamics. Therefore, it is widely believed that 
free-energy-driven dynamics, as it is relevant for protein 
folding, can reasonably be provided only by Boltzmann- 
Markov chains of conformational updates. One possibil- 
ity to achieve this is the application of a conventional 
Metropolis Monte Carlo method. Although the overall 
time scale is left open, this method allows for comparative 
folding (unfolding) studies at various temperatures. The 
kinetics of the folding (unfolding) transitions is thereby 
obtained by averaging over sufficiently many folding (un- 
folding) trajectories. 

In Fig. [Jj snapshots of a single GoL folding event 
of the heteropolymer with sequence S at T = 0.3 are 
shown at different times ^mc- The gray scale of the 
monomers encodes the variance of the monomer posi- 
tions, a^.. = x f — ^i 2 ; averaged over the MC time inter- 
val [£mc — 1 000, £mc + 1 000] . The first monomer is fixed 
and, therefore, not moving. The higher the mobility of 
a monomer in this time interval, the brighter is its color. 
For a%_. < 0.1, the monomer is rendered in black, and for 
<7^ > 1 hi white. The grayscales are linearly interpolated 
in-between these boundaries. Although there are periods 
of relaxation and local unfolding, a stable intermediate 
conformation is not present and the folding process is a 
relatively "smooth" process. This is also confirmed by 
the more quantitative analysis of the same folding event 
in Fig. [51 where the temporal averages of the similarity 
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FIG. 7: Single folding event for the heteropolymer with sequence S in the GoL model at T = 0.3 < Tf ~ 0.36. The description 
of the grayscale code is given in the text. 



parameters q and Q, and of the energy E are shown. 
Since the temperature lies sufficiently far below the fold- 
ing temperature (TV ~ 0.36), the free-energy landscape 
does not exhibit substantial barriers which hinder the 
folding process. 

Nonetheless, the chevron plot shown in Fig. [9] exhibits 
a rollover which means that the folding characteristics is 
not perfectly of two-state type, in which case the fold- 
ing (unfolding) branches would be almost linear [l6j |. 
In this plot, the temperature dependence of the mean 
first passage time t mfp is presented. We define t mfp 
as the average number of MC steps that are necessary 
to form at least 13 native contacts in the folding sim- 
ulations, starting from a random conformation. In the 
unfolding simulations, we start from the native state and 
tmfp is the number of MC steps required to reach a con- 
formation with less than 7 native contacts, i.e., 13 na- 
tive contacts are broken. In all simulations performed 
at different temperatures, tmfp is averaged from the 
first passage times of a few hundred respective folding 
and unfolding trajectories. Assuming a linear depen- 



dence at least in the transition state region, tmfp is di- 
rectly related to exponential folding and unfolding rates 
kf,u ~ 1/ t mfp ~ exp(— ef.u/ksT), respectively where 
the constants Sf >u determine the kinetic folding (unfold- 
ing) propensities. The dashed lines in Fig. [S] are tangents 
to the logarithmic folding and unfolding curves at the 
transition state. The slopes are the folding (unfolding) 
propensities and have in our case values of Sf w —1.32 
and e u ~ 5.0. 

In this variant of the chevron plot, which is similar to 
the presentations discussed in Refs. OHl], t ne tempera- 
ture T mimics the effect of the denaturant concentration 
that is in experimental studies the more generic exter- 
nal control parameter. The hypothetic intersection point 
of the folding and unfolding branches defines the transi- 
tion state. The transition state temperature estimated 
from this analysis coincides very nicely with the folding 
temperature Tf ~ 0.36 as identified in our discussion of 
the thermodynamic properties of the system. This re- 
sult also demonstrates that the description of the folding 
and unfolding transitions from the kinetic point of view 
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FIG. 8: Temporal averages of the native-contact and the over- 
lap similarity parameters q and Q, respectively, and the en- 
ergy for the folding event shown in Fig. [7] The temporal 
averages are calculated every 4000 MC steps over the time 
interval [t M c - 1 000, t M c + 1 000]. 
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FIG. 9: Chevron plot of the mean-first passage times from 
folding (•) and unfolding (o) events at different temperatures. 
The hypothetic intersection point corresponds to the transi- 
tion state. 



V. MESOSCOPIC HETEROPOLYMERS VS. 
REAL PROTEINS 

It is not the purpose of this manuscript to perform 
a coarse-grained analysis of the two-state folding char- 
acteristics of a specific, real protein. Rather, we have 
shown that - although also only exemplified for a sin- 
gle sequence - it is actually useful to study thermody- 
namic and kinetic properties of mesoscopic peptide mod- 
els without introducing atomic details. The main focus 
of such models is pointed towards general features of the 
folding transition (measured in terms of "order" parame- 
ters being specific for the corresponding transition, such 
as, e.g., the contact and overlap parameters investigated 
in our present study) that are common to a number of 
proteins behaving qualitatively similarly. It is then fur- 
thermore assumed that these proteins can be grouped 
into classes of certain folding characteristics. The se- 
quence S used in our paper is not obtained from a one- 
to-one hydrophobic-polar transcription of a real amino 
acid sequence. We do not think that such a mapping is 
particularly useful. Rather, S is considered as a repre- 
sentative that exhibits two-state folding characteristics 
in the coarse-grained model considered here. This im- 
plies, in general, that the classification of peptide fold- 
ing behaviors is not necessarily connected with detailed 
atomic correlations and particular contact-ordering. It 
is rather an intrinsic property of protein-like heteropoly- 
mers and can thus already be discovered employing mod- 
els on mesoscopic scales [19l [2fjj . 

Several proteins are known to be two-state folders and 
their folding transitions exhibit the features we have also 
seen in our present coarse-grained model study. A fa- 
mous example is Chymotrypsin Inhibitor 2 (CI2) one 
of the first proteins, where two-state folding character- 
istics has experimentally been identified. Clear signals 
of a first-order-like folding-unfolding transition were also 
seen in computational Go model analyses of that pep- 
tide [HI, [HI]. It is clear, however, that a precise character- 
ization of the transition state ensemble, which is required 
for a better understanding of the folding (or unfolding) 
process of a specific peptide (e.g., secondary-structure 
formation or disruption in CI2 [15|), is not possible. In 
the models used in our study, for example, only tertiary 
folding aspects based on hydrophobic-core formation are 
considered as being relevant. Nonetheless, as shown in 
our study of sequence S, macroscopic quantities or coop- 
erativity parameters manifest a qualitatively similar be- 
havior of the heteropolymer considered here, compared 
to real two-state folders. 



VI. SUMMARY 



is not only qualitatively, but even quantitatively consis- 
tent with the thermodynamic results from the canonical- 
ensemble analysis. 



In this study, we have analyzed the folding thermody- 
namics and kinetics of a knowledge-based hydrophobic- 
polar Go-like (GoL) model for heteropolymers at a meso- 
scopic scale. For a given sequence, the native conforma- 
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tion and the peak temperature of the specific heat as 
obtained from the physically motivated AB model were 
used to parametrize the GoL model. Since the chosen se- 
quence was known to exhibit two-state-like folding char- 
acteristics in the AB model, the definition of the cor- 
responding GoL model is appropriate for kinetic fold- 
ing studies, which are very time consuming with the AB 
model itself. The reason is that the Metropolis Monte 
Carlo method we used for the kinetic studies, is respon- 
sible for trapping effects and slowing down the Marko- 
vian dynamics of the physical model. Therefore, folding 
events are difficult to analyze in the AB model. This is 
the reason, why for statistical analyses of this model typ- 
ically general-ensemble methods are employed [HI, [28j . 
Unfortunately, these methods are not suitable to study 
kinetic aspects in a fixed-temperature ensemble. 

This work focuses on the question how kinetic as- 
pects of two-state folding behavior of realistic proteins 
can also be identified in a strongly simplified meso- 
scopic knowledge-based model, at least qualitatively. In 
Refs. [H, we could show that it is possible to reveal 
statistical folding properties of different folding charac- 
teristics employing the coarse-grained AB model. Here, 
we have found that the corresponding mesoscopic GoL 



model allows for the qualitative analysis of kinetic prop- 
erties of protein-like heteropolymers with two-state fold- 
ing behavior. Thermodynamic and kinetic properties 
were found to be even quantitatively consistent. 

The advantage of such simplified, coarse-grained mod- 
els is that they enable a more global, generalized view on 
the physics of conformational transitions accompanying 
protein folding processes. Our results show, in particular, 
that characteristic folding behaviors are not necessarily 
specific to microscopic details, but also an intrinsic prop- 
erty of hydrophobic-polar heteropolymers in general. 
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